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A  METHOD  FOR  ESTIMATING  THE  MECHANICAL  PROPERTIES 
OF  A  SOLID  MATERIAL  SUBJECTED  TO  INSONIFICATION  — 

PARTI:  THEORY 


INTRODUCTION 

Measurement  of  the  mechanical  properties  of  slab-shaped  materials  (i.e.,  plates)  is  important 
because  of  the  significant  contributions  made  by  these  properties  to  the  static  and  dynamic 
response  of  a  structure  under  various  loading  conditions. 

Resonant  techniques,1'3  used  to  identify  and  measure  Young’s  modulus  for  many  years,  have 
now  been  extended  to  measure  the  shear  modulus  of  materials.4  Such  procedures  are  based  on 
comparing  the  measured  eigenvalues  of  a  structure  to  those  predicted  from  a  model  of  the  same 
structure  formulated  with  well-defined  (typically  closed-form)  eigenvalues.  Resonant  techniques 
allow  measurements  at  resonant  frequencies  and,  in  general,  utilize  long,  thin  materials  so  that 
the  waves  propagate  only  in  a  longitudinal  or  torsional  direction.  The  methods  have  also  been 
extended  to  nonresonant  frequencies,  ’ '  which  typically  use  one  or  more  inverted  transfer 
functions  of  thin-bar  data  to  solve  for  Young’s  modulus.  Other  research  has  considered  the  bulk 
modulus  of  small-sized  materials  through  the  insonification  of  samples  and  the  use  of  lasers  to 
measure  the  resulting  strain.8 

Comparison  of  analytical  models  to  measured  frequency  response  functions  is  yet 
another  approach  used  to  estimate  the  stiffness  and  loss  parameters  of  a  structure.  When  the 
analytical  model  agrees  with  one  or  more  frequency  response  functions,  the  parameters  used  to 
calculate  the  analytical  model  are  considered  to  be  accurate.  However,  when  the  analytical 
model  is  formulated  with  a  numerical  method,  comparing  the  model  and  the  data  can  be 
difficult  because  of  the  dispersion  properties  of  the  materials.  Moreover,  if  a  single  type  of 
wave  motion  is  dominant  in  the  experiment,  the  parameters  associated  with  secondary  types  of 
wave  motion  are  unlikely  to  be  estimated  correctly. 

Some  investigators  have  concentrated  on  the  stiffness  and  loss  parameters  of  thin  plates,  with 
one  method  employing  transmission  zeros  to  calculate  unknown  coefficients  at  ultrasonic 
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frequencies  and  a  second  using  the  natural  frequencies  of  the  plate  to  determine  the  elastic 
constants.  In  addition,  plate  parameter  estimation  has  been  extended  to  measure  experimental 
plate  natural  frequencies  that  are  matched  to  a  numerical  simulation  for  estimating  elastic 
constants.15  Finally,  studies  of  plate  theory  that  relate  transmission  loss  and  reflection  loss  (echo 
reduction)  to  transfer  functions  across  the  plate  have  been  made,16  although  this  work  does  not 
contain  a  parameter  estimation  technique. 

A  theoretical  inverse  method  to  measure  complex,  frequency-dependent  dilatational  and 
shear  wavespeeds  that  are  present  in  thick,  isotropic  plates  is  developed  in  this  report.  After  the 
linear  equations  of  motion  of  the  system  are  derived  from  thick  plate  theory,  the  inverse  method 
is  used  to  combine  three  transfer  function  measurements  that  yield  closed-form  values  of  the 
complex,  frequency-dependent  Young  s  and  shear  moduli  and  complex,  frequency-dependent 
Poisson  s  ratio.  Simulations  are  presented  both  with  and  without  noise  to  show  how  the  method 
works  and  how  the  measurement  error  affects  the  analysis  results. 

A  typical  test  configuration,  as  seen  in  figure  1 ,  uses  a  speaker  to  project  acoustic  energy 
onto  the  material.  Intended  for  use  when  a  plate  is  to  undergo  acoustic  loading  on  its  face,  this 
technique  can  be  applied  to  automobile,  building,  and  submarine  materials. 


SYSTEM  MODEL 

As  described  above,  the  experiment  consists  of  insonifying  a  slab-shaped  test  specimen  by 
loading  the  structure  on  one  entire  side  with  an  acoustic  wave  originating  at  the  speaker.  The 
speaker  (or  projector)  is  located  at  such  a  distance  from  the  material  that  the  acoustic  wave 
nearly  becomes  a  plane  wave  by  the  time  it  has  established  contact  with  the  specimen.  This 
type  of  experiment  usually  conducted  at  multiple  frequencies  and  multiple  insonification 
angles  —  employs  a  frequency  sweep  (swept  sine)  at  three  different  insonification  angles  for  the 
inverse  method  presented  here.  During  insonification,  transfer  function  data  are  collected  from 
both  sides  of  the  plate,  either  with  accelerometers  that  measure  accelerations  or  with  laser 
velocimeters  that  use  reflected  light  to  measure  velocities.  In  the  swept-sine  mode,  transfer 
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Figure  1.  Test  Configuration  to  Insonify  and  Measure  Material 

functions  of  acceleration  divided  by  acceleration  or  velocity  divided  by  velocity  will  be  equal  to 
displacement  divided  by  displacement.  Finally,  the  time  domain  data  are  Fourier  transformed 
into  the  frequency  domain  and  then  recorded  as  complex  transfer  functions,  typically  with  a 
spectrum  analyzer. 

The  model  and  theoretical  analysis  are  based  on  the  following  assumptions: 

•  the  forcing  function  acting  on  the  plate  is  a  plane  wave  with  definite 
wavenumber  and  frequency  content; 
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•  the  corresponding  response  of  the  plate  is  at  a  definite  wavenumber  and 
frequency; 

•  motion  is  normal  and  tangential  to  the  plate  in  one  direction  (two- 
dimensional  system); 

•  the  plate  has  an  infinite  spatial  extent;  and 

•  the  particle  motion  is  linear. 

The  motion  of  the  test  specimen  shown  in  figure  1  is  governed  by  the  equation 


//V2u  +  (A  +  //)VV  •  u  =  P~^r  ,  (1) 

where  A  and  p  are  the  complex  Lame  constants  (  N/m2 ),  p  is  the  density  (  kg/m3 ),  t  is  time  (s), 

•  denotes  a  vector  dot  product,  and  u  is  the  Cartesian  coordinate  displacement  vector  of  the 
material. 

Figure  2  illustrates  the  coordinate  system  of  the  test  configuration.  Note  that  using  this 
orientation  results  in  b  =  0  and  a  having  a  value  that  is  less  than  zero,  with  the  thickness  of  the 
specimen,  h,  exhibiting  a  positive  value.  Equation  (1)  is  manipulated  by  writing  the 
displacement  vector  u  as 


ux(x,y,z,t) 
uy(x,y,z,t)  >  , 
u:(x,y,z,t) 


(2) 


where  x  is  the  location  along  the  plate  (m),y  is  the  location  into  the  plate  (m),  and  z  is  the 
location  normal  to  the  plate  (m),  as  shown  in  the  figure.  The  symbol  V  is  the  gradient  vector 
differential  operator  written  in  three-dimensional  Cartesian  coordinates  as 


xi  d  ■  d  ■  &  • 

v=;rl*+^+^- 

ox  dy  dz 


(3) 
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z  =  b  =  0 


Figure  2.  Coordinate  System  of  the  Test  Configuration 

with  ix  denoting  the  unit  vector  in  the  ^-direction,  i  denoting  the  unit  vector  in  the  ^-direction, 
and  i.  denoting  the  unit  vector  in  the  2-direction.  The  symbol  V2  is  the  three-dimensional 
Laplace  operator  applied  to  vector  u  as 


V2u  =  V2uxix  +  V2uyiy  +  V2mz4 


and  to  scalar  u  as 


-  d urvr  d uxv,  d 

V2w  =  V  •  Vw  = - ^  + ^  + - 22- 

'y'  'y'  dx2  dy2  dz1 


where  the  term  V  »u  is  the  divergence  and  is  equal  to 


_  du  du  dux 
V»u  =  — -  + — -  + — z- 
dx  dy  dz 


5 


The  displacement  vector  u  is  written  as 


u  =  V  <f>  +  V  x  (p  , 

where  ^  is  a  dilatational  scalar  potential,  x  denotes  a  vector  cross  product,  and  y/  is 
equivoluminal  vector  potential  expressed  as 


y/x(x,y,z,t) 

V  =  <  ysy(x,y,z,t)>  . 

y:(x,y,z,t) 


The  problem  is  formulated  as  a  two-dimensional  system,  and  thus  y  =  0,  uy(x,y,z,t)  =  0, 

and  d(-)/dy  =  0 .  Expanding  equation  (7)  and  breaking  the  displacement  vector  into  its 
individual  nonzero  terms  yields 


ux(x,z,t) 


d<f>(x,z,t) 


~.(v.  ,  <%(*■*>» 

dz  dx 


Equations  (9)  and  (10)  are  next  inserted  into  equation  (1),  which  results  ii 


|vV(x,z,/)  =  ^fe^ 

dt 2 


c)V2yr y{x,z,t) 


y(x,z,t) 
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where  equation  (11)  corresponds  to  the  dilatational  component  and  equation  (12)  corresponds  to 
the  shear  component  of  the  displacement  field.  Correspondingly,  the  constants  and  cs  are  the 
complex  dilatational  and  shear  wavespeeds,  respectively,  and  are  determined  by 


The  relationship  of  the  Lame  constants  to  the  Young’s  and  shear  moduli  is  shown  as 


(15) 


(16) 


where  E  is  the  complex  Young’s  modulus  (N/m2),  G  is  the  complex  shear  modulus  (N/m2),  and 
v  is  the  Poisson's  ratio  of  the  material  (dimensionless). 

The  conditions  of  infinite  length  and  steady-state  response  are  now  imposed,  allowing  the 
scalar  and  vector  potential  to  be  written  as 

(j){x,z,t)  =  <E>(z)exp(i£;c)exp(i<yO  (17) 

and 

y/y(x,z,t)  =  vP(z)exp(iA:x)exp(iryf) ,  (18) 
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where  i  is  the  square  root  of — 1,  eo  is  the  frequency  (rad/s),  and  k  is  the  wavenumber  with  respect 
to  the  x-axis  (rad/m).  Inserting  equation  (17)  into  equation  (11)  yields 


dz2 


j-^  +  or  O(z)  =  0  , 


where 


k2-k2 


k  -  — 

Kd  ~ 


Inserting  equation  (18)  into  equation  (12)  yields 


^  +  y8'4'(r)  =  0. 


where 


fi  =  Jk2-k2  , 


The  solution  to  equation  (19)  is 


O(z)  =  A(k,  o>)exp(i  az)  +  B(k, ry)exp(-iaz)  , 


and  the  solution  to  equation  (22)  is 
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'i'(z)  =  C(k,  a>)exp(ij3z)  +  D(k, co)exp(-ifiz)  , 


(26) 


where  A(k,co) ,  B(k,a>),  C(k,co) ,  and  D(k,co )  are  wave  response  coefficients  that  are 
determined  later  in  this  section.  Use  of  the  expressions  in  equations  (9)  and  (10)  now  allows  the 
displacements  to  be  written  as  functions  of  the  unknown  constants  as  follows: 

uz(x,z,t )  =  U,  (k,  z,  <w)exp(ifci')exp(i<yt) 

=  ^a\A{k,<a)Qxp(\az)  -  B(k,co)exp(-'mz  )]  (27) 

+  i£[C(A,<y)exp(i/?z)  +  D(k,a>)  exp(-i/fe)]}exp(iAx)exp(i®/) 
and 

ux{x,z,t)  =  Ux(k,z,(o)exp(ikx)exp(i(ot) 

=  §k\A(k,co)Qxp(iaz)  +  B(k,oj)exp{-\az)\  (28) 

-  i/?[C(&,  <x>)  exp(i/?z)  -  D(k,  a)  exp(-i/?z)]  }exp(i&x)  exp(irot)  . 

The  normal  stress  at  the  top  of  the  plate  ( z  =  b )  is  equal  to  the  opposite  of  the  pressure  load 
created  by  the  projector  that  acts  on  the  plate  surface  and  is  expressed  as 


= a + 2rt£^M+ .  -MxM) ; 

oz  ox 


the  tangential  stress  at  the  top  of  the  plate  is  zero  and  is  written  as 


=0. 

dz  dx 


The  normal  stress  at  the  bottom  of  the  plate  (z  =  a )  is  equal  to  zero,  yielding 
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(31) 


= a + 2 + 2^1^ = o 

dz  dx 


and  the  tangential  stress  at  the  bottom  of  the  plate  is  zero,  resulting  in 


The  applied  load  in  equation  (29)  is  an  acoustic  pressure  that  is  modeled  as  a  function  at 
definite  wavenumber  and  frequency  and  is  expressed  as 

Pv(x,z,t)  =  PQ((o)ex$(\kx)txv(\cof)  ,  (33 

where  the  wavenumber  k  is  found  using 

k  =  —  sin(0)  ,  nd 


in  which  cf  is  the  compressional  wavespeed  of  air  (m/s)  and  6  is  the  angle  of  incidence  of  the 
projector  with  the  2-axis  (rad). 

Assembling  equations  (l)-(34)  and  letting  b  =  0  yields  the  following  four-by-four  system  of 
linear  equations  that  models  the  structure: 

Ax  =  b  ,  (35) 

where  the  entries  of  equation  (35)  are 


An=-a2 X-2a2 //- Xk2  , 


■^12  ~  -^11  5 
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u> 

II 

N> 

(38) 

^14  —  ^13  ’ 

(39) 

A2l  =  -l/uka  , 

(40) 

A22  =  —A21  , 

(41) 

423=/'fl2  -M2  » 

(42) 

A  —A 

■^24  —  ^23  ’ 

(43) 

A31=Anexp(iaa)  , 

(44) 

A 32  =  v4nexp(-iaa)  , 

(45) 

^33  =  ~A13exp(ij3a)  , 

(46) 

A34  =  A13  exp(-i fid)  , 

(47) 

A41  =  A21  exp(i ad)  , 

(48) 

A42  =  -A21exp(-iaa)  , 

(49) 

A43  =  A23  cxp(i fid)  , 

(50) 

A44  =  A23exp(-i/?a)  , 

(51) 

II 

1 

/— s 

(52) 

62i  —  0  , 

(53) 

II 

o 

(54) 

o 

II 

(55) 
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From  equations  (35)-(55),  the  solution  to  the  constants  A{k,co),  B{k,co),  C(k,oj),  and  D{k,oS) 
can  be  calculated  at  each  specific  wavenumber  and  frequency  from 


Once  these  solutions  are  known,  the  transfer  function  between  the  wall  motion  in  the  z-direction 
atz  =  a  and  the  wall  motion  in  the  z-direction  at  z  =  b  is  written  in  closed- form  notation  from 
equations  (27)  and  (56)  as 

T(k,6))=  _ 4or/?A:2  sin(ar /;)  +  (/? 2  -  A:2)2  sm(/3h) _ 

U.(k,b,o))  4a^k2  sm(ah)cos(j3h)  +  (fi2  -  *2)2  cos(ah)sm(j3h)  '  ^ 

The  expression  in  equation  (57)  will  first  be  used  in  the  inverse  method  to  estimate  cc  and  then  to 
estimate  /?,  both  of  which  will  ultimately  be  employed  to  calculate  the  mechanical  properties  of 
the  material.  This  process  is  further  described  in  the  following  sections. 


INVERSE  MODEL  AT  ZERO  WAVENUMBER 

The  first  step  in  the  inverse  process  solves  for  the  response  at  zero  wavenumber  (typically 
referred  to  as  the  broadside  response)  to  determine  the  dilatational  wavespeed.  At  zero 
wavenumber,  the  angle  between  the  propagation  direction  of  the  insonification  energy  and  the 
z-axis  is  zero  and  the  response  of  the  structure  to  broadside  energy  is  comprised  entirely  of 
dilatational  waves  (i.e.,  no  shear  waves  are  excited  at  zero  wavenumber).  Moreover,  the  transfer 
function  given  in  equation  (57)  reduces  to 


7(0,  <y) 


1 

cos  (axh) 


(58) 


where  7,  (or  Rx)  represents  the  data  from  the  experiment  with  an  insonification  angle  of  zero 
(typically  a  frequency-dependent,  complex  number)  and  the  subscript  1  denotes  the  first 
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experiment.  Equation  (58)  can  be  expanded  into  real  and  imaginary  parts  and  solved,  resulting 
in  a  value  for  al  at  every  frequency  for  which  a  measurement  is  made.  The  solution  to  the  real 

part  of  is 


— Arccos(s)  +  — 
2  h  2  h 


Re(«1)  = 


— Arc  cos(-^)  +  — 
^2  h  2  h 


n  even 


n  odd 


(59) 


where 

s  =  \R e(^)]2  +[Im(i?1)]2 -V{[R e(^)]2  +[Im(/?1)]2f  -{2[Re(i?1)]2 -2[Im(^)]2  -l}  , 

(60) 

n  is  a  non-negative  integer,  and  capital  A  denotes  the  principal  value  of  the  inverse  cosine 
function.  The  value  of  n  is  determined  from  the  function  s,  which  is  a  periodically  varying 
cosine  function  with  respect  to  frequency.  At  zero  frequency,  n  is  0.  Every  time  s  cycles 
through  n  radians  (180°),  n  is  increased  by  1 .  When  the  solution  to  the  real  part  of  al  is  found, 

the  solution  to  the  imaginary  part  of  ax  is  written  as 


Im(«i)  =  |loge 
h 


RgCgi) 

cos[Re(ttj  )h] 


ittW  | 

sinlRe^j)^]] 


(61) 


The  real  and  imaginary  parts  of  ax  from  equations  (59)  and  (61),  respectively,  are  combined  to 
yield  the  complex  wavenumber.  Because  this  measurement  is  made  at  zero  wavenumber 
(k  =  0),  a{  is  equal  to  the  dilatational  wavenumber.  Thus,  the  dilatational  wavespeed  is 
expressed  as 


co 

p.e(a1)  +  ilm(a1)] 


(62) 
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To  solve  for  the  shear  wavespeed,  the  specimen  must  be  excited  at  a  nonzero  wavenumber  as 
shown  in  the  next  section. 


INVERSE  MODEL  AT  NONZERO  WAVENUMBER 

The  next  step  solves  for  the  response  at  nonzero  wavenumber  to  determine  the  shear 
wavespeed.  The  transfer  function  at  the  nonzero  wavenumber  is  given  in  equation  (57)  and  can 
be  expressed  for  the  nonzero  angle  of  insonification  as 


T(k,co)  = 


_ 4Q'2/?2ft22  sin(o;2ft)  +  -  *22)2  sin(/?2/;) _ 

4  a2P2k\  sin(or2/?)cos(/?2/i)  +  (, 02  -  k\f  cos(a2h)sm(J32h) 


-  T  = 
12 


where  T2  (or  R2 )  represents  the  data  from  the  experiment  at  the  nonzero  insonification  angle 
(typically  a  frequency-dependent,  complex  number)  and  the  subscript  2  denotes  the  second 
experiment.  It  is  noted  that  a2  in  equation  (63)  is  different  from  a,  in  equation  (58),  with  this 
difference  based  on  the  k 2  term  shown  in  equation  (20),  where  the  wavenumber  a  is  defined. 
Due  to  the  complexity  of  equation  (63),  there  is  no  simple  method  for  rewriting  the  equation  as 
a  function  of  / ?2 ,  which  is  the  variable  that  is  to  be  estimated.  Rather,  equation  (63)  is 
rewritten  as 

/(&)  =  0  =  4 a2p2k]  sm(a2h) [cos(/?2 h) -  R2]  +  (J32  -  k2)2  s\n(p2h)[cos{a2h)  -  R2]  , 

(64) 

where  the  problem  now  becomes  finding  the  zeros  of  the  right-hand  side  of  equation  (64)  or,  in 
the  presence  of  actual  data  that  contain  noise,  finding  the  relative  minima  of  the  right-hand  side 
of  this  equation,  with  the  purpose  of  determining  which  relative  minimum  corresponds  to  shear 
wave  propagation  and  which  relative  minima  are  extraneous.  Because  equation  (64)  has  a 
number  of  relative  minima,  zero-finding  algorithms  are  not  applied  to  this  function,  as  they 
typically  do  not  fmd  all  the  minima  locations  and  are  highly  dependent  on  initial  guesses.  The 
best  method  for  locating  all  the  minima  is  to  plot  the  absolute  value  of  the  right-hand  side  of 
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equation  (64)  as  a  surface,  with  the  real  part  of  /?2  on  one  axis  and  the  imaginary  part  of  /?2  on 
the  other  axis.  The  value  a2  is  determined  from 

a2  =  A2  ~kl  =  V«i2  ~kl  >  (65) 

so  that  equation  (64)  is  a  function  of  only  (complex)  /?2 .  Once  this  function  is  plotted,  the 
minima  can  be  easily  identified  and  the  corresponding  value  of  (/?2)m  at  the  location  of  the 

minima  can  be  determined  by  examination  of  the  minimum  location  —  an  approach  that  is 
sometimes  referred  to  as  the  grid  method.  The  shear  wavenumber(s)  and  wavespeed(s)  are  then 
determined  from 


(6«) 


and 


(c,)m 


0) 

(*X’ 


(67) 


where  the  subscript  m  denotes  each  minima  value  that  was  found  by  inspecting  the  surface 
formed  from  equation  (64).  Determination  of  the  correct  index  of  m  that  corresponds  to  shear 
wave  propagation  is  presented  in  the  next  section. 


DETERMINATION  OF  PROPERTIES  FROM  WAVESPEEDS 


The  process  of  determining  the  material  properties  from  the  wavespeeds  begins  with 
calculation  of  the  Lame  constants  from  equations  (13)  and  (14),  which  is  expressed  as 


Pm  =  P(csfm 


(68) 


and 


4  =  Pc2d~2p(cs)2m  . 


(69) 
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To  determine  the  correct  index  m  that  corresponds  to  the  actual  wave  propagation  rather  than  to 
an  extraneous  solution,  a  third  set  of  measurements  is  made  at  a  nonzero  incidence  angle  that  is 
not  equal  to  the  angle  used  in  the  previous  section.  That  is,  the  model  in  equation  (63)  is 
calculated  from  the  estimated  material  properties,  with  a  residual  value  established  from  this 
third  set  of  measurements.  Each  m-indexed  residual  at  a  specific  frequency  is  defined  as 


(*3)„  = 


_ 4a3 (A )m *3  sm(a3h)  +  [(/?3 )2m-klf  sin[(/?3 )m  h] _ 1_ 

4a3  )m  *3  sm(cc3h) cos[(/?3 )m  h]  +  [(/?3  fm  -  kl ]2  cos (a3h) sin[(/?3  )m  h\  R3  ’ 


where 


^3  —  \  k^  k->  —  —  k? 


(A)m  -  V(A)m  +  *2  ~k]  ,  ^72) 

with  the  subscript  3  denoting  the  third  experiment.  The  smallest  residual  value  corresponds  to 

the  correct  value  of  index  m  and  the  correct  values  of  the  Lame  constants.  Poisson’s  ratio  is  then 
calculated  using 


2(//  +  A) 


Young’s  modulus  can  be  calculated  from 


£_  2//(2//  +  3A) 

2(M  +  A)  ’ 


and  the  shear  modulus  from 
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NUMERICAL  SIMULATION 


The  inverse  measurement  method  presented  here  can  be  simulated  by  means  of  a  numerical 
example  that  uses  a  soft  rubberlike  material  with  the  following  properties:  Young’s  modulus  E  of 
[(le8  -  i2e7)  +  (5e3 /  -  i3e2/)]  N/m2 ,  where /  is  frequency  in  Hertz;  Poisson’s  ratio  v  equal  to 
0.40  (dimensionless);  density  p  equal  to  1200  kg/m3;  and  thickness  h  of  0.1  m.  A  compressional 
(acoustic)  wave  velocity  of  cf  =  343  m/s  for  air  was  used.  All  other  parameters  can  be  calculated 
from  these  values. 

Insonification  angles  of  0°,  20°,  and  40°  are  chosen  to  illustrate  the  technique,  with  three 
different  amounts  of  additive  noise  used  during  the  simulation:  no  additive  noise,  1%  additive 
noise,  and  2%  additive  noise.  Additive  noise  is  included  in  the  transfer  function  using  the 
equation 

Te  (co)  =  T(co)  +  e{R e[T(a)]aa  +  ilm[T(co)]ab }  ,  (76) 


where  e  is  the  amount  of  additive  noise  added  to  the  transfer  function  and  <7a  and  crb  are  random 
numbers  with  zero  mean  and  a  variance  of  one.  The  value  e  is  also  called  the  transfer  function 
error  value,  and  it  represents  the  deviation  of  the  transfer  function  from  perfect  data. 

Figure  3  is  a  plot  of  the  transfer  functions  of  equation  (57)  at  0°,  20°,  and  40°  insonification 
angles  versus  frequency.  The  top  plot  illustrates  the  magnitude  and  the  bottom  plot  presents  the 
phase  angle.  Although  the  transfer  functions  shown  are  without  additive  noise,  the  transfer 
functions  with  additive  noise  were  similar  to  those  seen  in  this  figure.  Once  the  transfer 
functions  are  known  (typically  by  measurement,  but,  in  this  report,  by  numerical  simulation),  the 
dilatational  wavespeed  can  be  estimated  from  equations  (59)-(62). 

Figure  4  is  a  plot  of  the  function  s  versus  frequency  without  additive  noise,  and  figure  5 
illustrates  the  actual  dilatational  wavespeed  and  the  estimated  dilatational  wavespeed  without 
additive  noise,  with  1%  additive  noise,  and  with  2%  additive  noise  versus  frequency.  The 
top  plot  is  the  real  component  and  the  bottom  plot  is  the  imaginary  component  of  the 
wavespeed. 
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Figure  6  plots  the  absolute  value  of  the  magnitude  of  the  surface  defined  in  equation  (64) 
versus  the  real  and  imaginary  components  of  /?2  at  1800  Hz  without  additive  noise.  The  surface, 
discretized  with  400  by  400  points,  has  its  magnitude  expressed  in  a  decibel  scale.  Seven  distinct 
local  minima  exist,  with  the  first  local  minimum  corresponding  to  /?2  =  0,  which  is  a  value  that 
implies  no  shear  wave  propagation  (a  physically  unrealizable  condition  at  nonzero  wavenumber). 
The  other  six  local  minima  are  labeled  on  the  figure  with  bold  numbers  and  are  processed  at  a 
third  measurement  location  according  to  equation  (70). 

The  results  for  no  additive  noise  are  provided  in  table  1 .  As  is  evident,  local  minimum 
number  3  has  the  smallest  residual  value  and  corresponds  to  the  shear  wave  propagation.  The 
value  for  (/?2)3  is  equal  to  60.9  +  5.9i  compared  to  the  actual  value  of  /?2  which  is  61.0  +  5.9i. 
The  small  difference  between  these  two  values  can  be  attributed  to  the  discretization  of  the 
surface  as  shown  in  the  figure.  Additionally,  it  is  noted  that  the  local  minima  numbered  4,  5,  and 
6  have  an  imaginary  value  almost  equal  to  zero,  which  would  correspond  to  a  loss  factor  of  zero. 
Most  soft  materials  have  moderate  to  large  loss  values,  which  means  that  these  minima  could  be 
disregarded  based  on  this  criterion.  The  variable  p2  at  1800  Hz  was  next  estimated  using  1% 
and  2%  additive  noise,  as  shown  in  tables  2  and  3,  respectively.  These  results  illustrate  that  the 
residual  method  is  effective  when  noise  is  present  in  the  data. 

Figure  7  is  a  plot  of  the  actual  shear  wavespeed  and  the  estimated  shear  wavespeed  without 
additive  noise,  with  1%  additive  noise,  and  with  2%  additive  noise  versus  frequency.  The  top 
plot  presents  the  real  component  and  the  bottom  part  shows  the  imaginary  component  of  the 
wavespeed.  The  figure  was  constructed  by  repeating  the  measurement  process  shown  in  figure  6 
at  each  specified  frequency. 

Figure  8  plots  the  actual  shear  modulus  and  the  estimated  shear  modulus  without  additive 
noise,  with  1%  additive  noise,  and  with  2%  additive  noise  versus  frequency.  The  top  plot 
represents  the  real  component  and  the  bottom  part  shows  the  loss  factor  of  the  shear  modulus. 

Figure  9  is  a  plot  of  the  actual  Young’s  modulus  and  the  estimated  Young’s  modulus 
without  additive  noise,  with  1%  additive  noise,  and  with  2%  additive  noise  versus  frequency. 
The  top  plot  is  the  real  component  and  the  bottom  part  is  the  loss  factor  of  the  Young’s 
modulus. 
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Figure  6.  Magnitude  of  Surface  Versus  Real  and  Imaginary  /?2  at  1800  Hz 


Table  1.  Values  of  (J32)m  and  (c3)„,  at  the  Local  Minima 
Without  Additive  Noise 


Local  Minima 

Number  m 

Value  of 

(A). 

Residual  (f3)m 
(Equation  (70)) 

1 

22.1  +  5.2i 

0.302 

2 

38.4  +  2.5i 

2.006 

3 

60.9  +  5.9i 

0.001 

4 

94.3  +  0.5i 

0.506 

5 

125.1  +0.6i 

0.335 

6 

157.4  +  0.5i 

0.378 
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Table  2.  Values  of  (/?2)m  and  {s3)m  at  the  Local  Minima 
with  1  %  Additive  Noise 


Local  Minima 
Number  m 


1 


2 


3 


Value  of 

(A)m 


22.1  +  5.2i 


39.3  +  2.2i 


60.9  +  6.3i 


94.3  +  0.5i 


Residual  (e3)m 
(Equation  (70)) 


0.275 


2.170 


0.045 


0.516 


5 

125.0  +  0.6i 

0.312 

6 

157.4  +  0.5i 

0.371 

Table  3.  Values  of  (p2)m  and  (e3)m  at  the  Local  Minima 
with  2%  Additive  Noise 


Local  Minima 
Number  m 


Value  of 

(A). 


22.1  +  5.3i 


38.4  +  2.7i 


60.9  +  5.2i 


94.3  +  0.5i 


125.0 +  0.5i 


157.4 +  0.5i 


Residual  {s3)m 
(Equation  (70)) 


0.344 


1.979 


0.009 


0.519 


0.361 


0.398 
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Figure  8.  Actual  and  Estimated  Shear  Modulus  Versus  Frequency 
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Figure  9.  Actual  and  Estimated  Young’s  Modulus  Versus  Frequency 


Figure  10  illustrates  the  actual  real  Poisson’s  ratio  and  the  estimated  real  Poisson’s  ratio 
without  additive  noise,  with  1%  additive  noise,  and  with  2%  additive  noise  versus 
frequency.  Because  the  numerical  example  is  formulated  using  a  Poisson’s  ratio  that  is 
strictly  real,  no  imaginary  component  is  shown  in  this  plot.  Imaginary  values  of  Poisson’s 
ratio  are  possible  and  have  been  shown  to  theoretically  exist.17 

Finally,  the  additive  noise  versus  the  parameter  estimation  error  is  shown  in  table  4. 
Parameter  estimation  error  was  computed  from 


e=  1  yK„K)~^K)l 


where  6  is  the  parameter  estimation  error,  tcest  (<y„)  is  the  estimated  value  of  the  parameter  at 
the  nth  frequency  value,  Kact  (con  )  is  the  actual  value  of  the  parameter  at  the  nth  frequency 
value,  and  N  is  the  total  number  of  frequencies  at  which  an  estimate  was  computed.  It  is  noted 
that  at  very  low  frequencies  the  routine  did  not  produce  a  realistic  estimate  of  the  parameters, 
and  thus  these  estimates  are  not  included  in  the  table. 

Table  4  shows  that  all  the  parameters  are  being  estimated  accurately  with  the  inverse 
technique.  At  1%  additive  noise,  the  maximum  error  in  the  wavespeeds  is  2.6%  and  the 
maximum  error  of  the  material  properties  is  5.1%.  At  2%  additive  noise,  the  maximum  error  in 
the  wavespeeds  is  3.4%  and  the  maximum  error  of  the  material  properties  is  6.8%.  Poisson’s 
ratio  has  a  parameter  estimation  error  of  1.4%  when  the  transfer  functions  have  2%  additive 
noise,  which  is  remarkably  low  for  a  parameter  that  is  historically  very  difficult  to  measure. 
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Figure  10.  Actual  and  Estimated  Poisson’s  Ratio  Versus  Frequency 


CONCLUSIONS 


In  this  report,  a  theoretical  method  for  estimating  the  mechanical  properties  of  slab-shaped 
materials  subjected  to  insonification  has  been  derived.  Measurement  of  the  transfer  functions  of 
acceleration  across  the  plate  at  three  different  insonification  angles  shows  that  the  dilatational 
wavespeed,  shear  wavespeed,  Young’s  modulus,  shear  modulus,  and  Poisson’s  ratio  can  be 
accurately  determined.  Moreover,  the  technique  is  relatively  immune  to  noise  mechanisms  that 
are  sometimes  present  in  the  measurements. 

It  is  recommended  that  actual  measurement  data  be  used  to  evaluate  the  effectiveness  of  the 
inverse  method. 
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